Liverworts show a globally consistent mid‐elevation richness peak

Abstract The study of elevational gradients allows to draw conclusions on the factors and mechanisms determining patterns in species richness distribution. Several earlier studies investigated liverwort diversity on single or few elevational transects. However, a comprehensive survey of the elevational distribution patterns of liverwort richness and their underlying factors is lacking so far. This study's purpose was to fill this gap by compiling an extensive data set of liverwort elevational patterns encompassing a broad diversity of mountains and mountain ranges around the world. Using polynomial regression analyses, we found a prevalence of hump‐shaped richness patterns (19 of 25 gradients), where liverwort species richness peaked at mid‐elevation and decreased towards both ends of the gradient. Against our expectation and unlike in other plant groups, in liverworts, this pattern also applies to elevational gradients at mid‐latitudes in temperate climates. Indeed, relative elevation, calculated as the percentage of the elevational range potentially inhabited by liverworts, was the most powerful predictor for the distribution of liverwort species richness. We conclude from these results that the admixture of low‐ and high‐elevation liverwort floras, in combination with steep ecological gradients, leads to a mid‐elevation floristic turnover shaping elevational patterns of liverwort diversity. Our analyses further detected significant effects of climatic variables (temperature of the warmest month, potential evapotranspiration, and precipitation of the warmest month) in explaining elevational liverwort richness patterns. This indicates that montane liverwort diversity is restricted by high temperatures and subsequent low water availability especially towards lower elevations, which presumably will lead to serious effects by temperature shifts associated with global warming.


| INTRODUC TI ON
Biodiversity is not equally distributed in space, neither latitudinally nor elevationally. The latitudinal gradient in species richness was the first ecological pattern identified and was already described by Alexander von Humboldt in the early 19th century (Hawkins, 2001).
The decrease in species richness from low latitudes near the equator towards higher ones at the poles has been found to be a general pattern applying to most major taxonomic groups (e.g., foraminifera (Culver & Buzas, 2000), bats and marsupials (Lyons & Willig, 2002), vascular plants (Mutke & Barthlott, 2005)), including liverworts (Wang et al., 2017).
Quite independently of latitude, mountains are global hotspots of biodiversity (Körner, 2004). Their steep climatic gradients, together with the heterogeneity of terrain and topographical structure within short distances lead to a great variety of habitats. These promote several processes associated with speciation and maintenance of species diversity (Perrigo et al., 2020). Mountain biodiversity shows a diverse range of elevational richness patterns. In parallel to the latitudinal diversity gradient, it was long believed that species richness would mostly decline linearly with increasing elevation, until Rahbek (1995) revealed that on the majority of elevational gradients richness peaks at some intermediate point of the gradient, although other patterns (linear decreases; constant and then declining; and even increases) also occur (McCain & Grytnes, 2010;Rahbek, 2005). Nevertheless, the hump-shaped pattern has been found to be prevalent in most organisms (Rahbek, 2005). For instance, Grytnes (2003b) found a hump-shaped relationship between vascular plant species richness and elevation along five of seven study transects in Norway, and declining richness on the other two.
For ferns on a global scale,  also found the majority of their 20 study transects to show a hump-shaped richness pattern for ferns, but increases, decreases, and constant richness were also found.
Of course, elevation is not by itself the determining factor of species richness (Gaston, 2000). Rather, it is an indicator of a mélange of biotic and abiotic determinants, which influence species richness. Determinants to be considered are manifold (Vetaas, 2021) and include variables relevant for metabolism like the availability of energy and water such as temperature, and humidity (Kluge et al., 2017) and other environmental factors like soil properties (Liu et al., 2020;Ohdo & Takahashi, 2020;Sánchez-González & López-Mata, 2005), land surface area (Karger et al., 2011), past geological processes including orogeny and plate tectonics (Descombes et al., 2017;Hagen et al., 2019;Zhao & Li, 2017), phylogenetic processes such as diversification rates and phylogenetic diversity Scholl & Wiens, 2016), biogeographic processes including source-sink effects and overlaps of lowland and montane floras (Grytnes, 2003a;Kessler, Hofmann, et al., 2011). One of the major challenges in identifying these drivers of elevational richness patterns is that many potential explanatory factors covary with elevation and among each other. This is especially relevant when species richness is studied along a single transect (Kluge et al., 2006).
Thus, the comparative analysis of multiple transects, which preferably exhibit different parameter peculiarities such as different climatic conditions and elevational extents, offers the opportunity to bypass this issue of covariation Lomolino, 2001;McCain, 2009;Vetaas et al., 2018).
Liverworts have evolved as one of the most basal groups of land plants, and comprise between 5000 and 7500 extant species (Gradstein, 2013;Von Konrat et al., 2010). Mostly being small to tiny plants, liverworts have evolved adaptations where their rhizoids allow them to grow adherent on impervious and steep surfaces, enabling them to inhabit a large variety of different habitats in addition to soil (Proctor et al., 2007) such as rocks, dead wood, tree trunks, branches, and even leaf surfaces of vascular plants (Jones & Dolan, 2012). Due to the lack of thick cuticles, liverworts desiccate quickly in the absence of liquid water and under low humidity, while suspending their metabolism. The ability to recover from desiccation (poikilohydry) is species-specific, and its success depends on multiple conditions (e.g., duration of desiccation and humidity; León-Vargas et al., 2006). Generally speaking, liverwort species inhabiting exposed and dry surfaces are most adapted to desiccation, followed by epiphytes of branches and the canopy. Species inhabiting tree bases and forest floor are less tolerant to drying (Franks & Bergstrom, 2000;Proctor et al., 2007), and aquatic bryophyte species are the least drought tolerant (Glime, 2017). In general, cool and humid conditions are most suitable for liverworts. Such conditions are often found at intermediate elevations of mountains (Körner et al., 2011;Rahbek, 1995).
There are various studies providing evidence for a hump-shaped elevational richness gradient in liverworts from different parts of the world, most of them considering individual transects (Henriques et al., 2016;Tabua et al., 2017;Wolf, 1993), although some are on a regional scale (e.g., Grau et al., 2007), or compare two local data sets (Ah-Peng et al., 2012;Iskandar et al., 2020). But divergent patterns of the elevation-dependent change in liverwort species richness have also been described. For instance, Bruun et al. (2006) reported increasing diversity in Fennoscandia. Liverwort diversity in general has been linked to both local-scale predictors including microhabitat variation (Ah-Peng et al., 2007), disturbance effects (Boch et al., 2019), and bark roughness (Gradstein & Culmsee, 2010), as well as to regional-scale factors, in particular humidity and temperature, and fog occurrence (Aranda et al., 2014;Frahm & Ohlemüller, 2001;Gehrig-Downie et al., 2013;Henriques et al., 2016). However, a comprehensive comparative analysis of elevational gradients in species richness of liverworts is still lacking, so that no general assessment is available to date.
In this study we intended to close this knowledge gap by studying liverwort species richness on a broad diversity of mountains and mountain ranges at a global scale to address the following questions: 1. What kind of elevational richness patterns are most common for liverworts? 2. Which parameters can best explain these overall elevational richness patterns of liverworts when analyzed in combination? 2 | MATERIAL S AND ME THODS

| Data sources
We pooled the liverwort species richness per elevation of 50 data sets from 31 published projects/expeditions and one unpublished transect studied by us (Table 1; Appendix S1). Data obtained via various study methods were included, provided that they comprise a minimum of four elevational steps (or could be transformed into such). Furthermore, the original data must contain at least 10 taxa in total, to allow for the development of an elevational pattern in species richness. We avoided data sets in which the land surface area of mountains did not successively decline with elevation (Körner, 2000), to avoid atypical area effects caused, e.g., by high plateaus. We analyzed the species richness of 25 gradients totaling 445 taxa per elevation records from all over the world for liverwort species richness in general (i.e., without specified habitat). In addition, we conducted analyses from 16 and nine locations (163 and 90 records) for species found to inhabit epiphytic and non-epiphytic habitats, respectively ( Figure 1). Our aim was to analyze the liverwort communities at a regional rather than local scale and data on the elevational distribution of species were sometimes taken from broader, regional works. For example, elevational data on species recorded from Napo province, Ecuador, were extracted from the Ecuadorian checklist (Gradstein, 2020), even though the information in the source publication referred to the wider northern Andean region and was not restricted to the Napo gradient. The Ecuadorian data were subsequently corrected for the species with lower elevational records elsewhere than in Napo alone based on personal communication with S.R. Gradstein.

| Explanatory variables: elevation and climate
In order to study the relationships between species richness and elevation, as well as different climatic explanatory variables, we used the relative species richness (relRich) of each location as a response variable to give equal weight to all transects and to minimize the influence of sampling intensity and area of the individual transects. relRich was calculated by dividing the species richness of each plot or elevational belt of a given transect by the maximum species richness per plot/belt of that transect. We then transformed relRich according to Smithson and Verkuilen (2006) to fit the requirements of beta regression models.
Elevation was included both as absolute elevation (absElev) and as relative elevation (relElev) in our analyses. absElev provides information on the environmental conditions that change in relation to elevation, whereas relElev places the liverworts richness in relation to the elevational extent of the transect, which can influence richness via source-sink effects (Kessler, Hofmann, et al., 2011), overlaps between lowland and highland floras (Grytnes, 2003a), and island effects on small mountain ranges (Kessler & Kluge, 2022).
For this reason, relElev no longer reflects direct relationships with climate as does absElev but rather informs about independent biogeographical processes, which in contrast to climatic factors would also explain hump-shaped elevational richness patterns in liverworts on mountains with short elevational extent. The quadratic term (relElev 2 ) would then be expected to be a highly significant predictor for species richness. We calculated relElev of each transect by defining the elevational range potentially inhabited by liverworts, from sea level to the highest summit of each transect, unless there was a permanent snow line at highest elevations, in which case this was taken as the uppermost limit. relElev was then calculated as the percentage of the upper boundary of the respective gradient.
In the absence of exact GPS coordinates for most data points of published studies, we obtained the climate data by determining 7-12 likely coordinates in each mountain range and the respective elevations per gradient from Google Earth. These coordinates were located in the respective region of interest, as close to the original plots as it was assessable from the source publications. We then obtained Bioclim data for temperature (annual mean (Bio1), max of warmest month (Bio5), min of coldest month (Bio6)) and precipitation amount (annual (Bio12), of wettest month (Bio13), of driest month (Bio14), of warmest quarter (Bio18), and of coldest quarter (Bio19)), as well as potential evapotranspiration (PET) from CHELSA (Karger et al., 2017) for these coordinates. We selected these factors because mean parameters inform about the average conditions at a given elevation, whereas the maximum and minimum parameters approximate the likely biologically limiting factors, e.g., when the ranges of some species are determined by frost. PET was included as the maximum PET value of the year (maxPET) in our analyses. Because temperature decreases with increasing elevation (Spellman, 2013), we predicted temperatures and maxPET for our data set using linear models. The relationship between precipitation and elevation is more complex (Karger et al., 2017;McCain & Grytnes, 2010). Hence, we fitted linear or polynomial models (quadratic or cubic) for the precipitation data. We only implemented predicted precipitation variables from models with (adjusted) R 2 > .5 in our data set. Two study regions (i.e., Switzerland and Nepal) were excluded a priori because we assumed precipitation amounts to be too heterogeneous to predict with our method. For locations where detailed GPS coordinates were available (Uganda, Madagascar, Mt. Wilhelm on Papua New Guinea, Madeira), we included the according precipitation data directly, without employing models.
For the interpretation of causality, we used absElev to describe the elevational richness patterns but without assuming any direct causal influence of elevation. Rather, we used climatic factors to infer the influence of physiological limitations of the liverworts in determining their richness patterns. Finally, we used relElev as an indicator of biogeographical processes related to mountain height, independently of climate.

TA B L E 1 Properties of source data included in this study.
No.

| Data analyses: best predictor combinations
In the first step of analyses, we related the distribution of relRich to elevation and the climatic factors using beta regression models, allowing for both linear and unimodal relationships. To assess the goodness of model fit we used the pseudo-R 2 (R 2 p ) function suggested by Ferrari and Cribari-Neto (2004) implemented in the R package betareg 3.1.-4 (Cribari-Neto & Zeileis, 2010). Likelihood ratio tests of nested models were conducted to decide if first or second order regression models were preferable.
In order to explain the highest possible variance of global elevational richness distribution in liverworts, we combined elevational and temperature-related predictors with pairwise Pearson's correlation coefficients <.7 after centering the predictors (i.e., subtracting the mean). We also included quadratic terms of all explanatory variables (Model1). To test the influence of precipitation on the largest possible subset of species richness, we included the precipitation variable with the best fit in the simple regression (Model2). For the sake of comparability, we additionally regressed the reduced species richness data used for Model2 against combinations of elevational and temperature-related predictors without precipitation (Model3).
The combined full models were ranked according to the corrected Akaike information criterion (AICc; Hurvich & Tsai, 1989), and are referred to as "best model" within ΔAICc = 2. The goodness of model fit was again assessed using R 2 p . Subsequently, we reduced the predictors of the best models stepwise until all predictors were significant (p < .001) to avoid overfitting. In order to unravel the effects from the individual versus the joint effects of predictor combinations, we conducted variance partition analyses based on partial regression: R 2 p full model = R 2 p predictor 1 + R 2 p predictor 2 + R 2 p shared .
All computational analyses were conducted in R 4.0.2 (R Core Team, 2017) with betareg 3.1.-4 (Cribari-Neto & Zeileis, 2010) to fit the beta regression models, lmtest 0.9-38 (Zeileis & Hothorn, 2002) to conduct likelihood ratio tests of nested models, VennDiagram 1.6.0 (Chen, 2018) to draw the Venn diagrams, and the raster package 3.4-5 (Hijmans, 2020) for the extraction of the Bioclim data.  Three locations (Itatiaia NP in Brazil, Cordillera Oriental in Colombia, and Peruvian Andes) did not show a significant elevational richness pattern. The analysis of the combined regional data resulted in a unimodal relationship between relRich and absElev with an R 2 p value of 0.19 (p < .001; Figure 2; Table 2).

| Simple regressions
Global relRich also showed a unimodal relationship to relElev, with a maximum at almost exactly the center of the global gradient showed a weak hump-shaped relationship between relRich and Bio14 (p lrtest <.1). Bio14 was not significant on the remaining elevational gradients. On a global scale, species richness increased linearly with increasing precipitation, which accounted for 17% of the variance (p < .001; Figure 2; Table 2).
The stepwise reduction of the best models resulted in a model with slightly lower R 2 p (.39), which comprised three predictors: relElev 2 and maxPET + maxPET 2 alone accounted for 23%, and 13% (Figure 3), respectively. The joint effect of the two variable fractions explained an additional 3%.
Regressing the reduced species richness data set against all possible combinations of relative elevation, and temperature-related variables plus precipitation of the driest month (Bio14), resulted in 10 best models within ΔAICc ≤2 (Model2). relElev, relElev 2 , Bio5 2 and Bio14 were constant components (p < .01) of all models. MaxPET 2 and Bio6 were present in nine and six best models (Appendix S4).
Further components of the best models were maxPET (five models), and Bio1 2 and Bio6 2 (three models each). The mean R 2 p -value of the best models was .48. After stepwise reduction of the predictors, the resulting model comprised relElev, relElev 2 , and Bio14. The R 2 p -value of the model slightly declined to .45, of which both terms of relative elevation accounted for the largest proportion (27%), and Bio14 for 17%. The joint effect of both fractions was negligible (1%).
Combining elevation variables with climatic factors resulted in 12 and seven (including Bio14) best models with R 2 p -values up to .21 (epiphytes), and seven and five (including Bio14) best models and R 2 p -values up to .21 (non-epiphytes). Stepwise reduction of the predictors was not possible to the significance level of p < .001 for the remaining predictors without reducing R 2 p -values to a fractional amount of the more complex model, with one exception: the reduced data set of non-epiphytic relRich (Models 2 and 3) yielded in a model containing relElev 2 as the only predictor (p < .001) with a noteworthy R 2 p value of .18.

Precipitation of driest month [mm]
et al., 2020; Tusiime et al., 2007). Our study confirms that, overall, the hump-shaped pattern dominates, and that the richness peak is on average found in the middle of the elevational gradients. The prevalence of the hump-shaped pattern in epiphytic species richness (56%), in contrast to non-epiphytes (33%), suggests that the pattern in overall species richness is largely driven by the epiphytes.
In this study, we used a set of potential predictors to explain these elevational richness patterns. In accordance with many previous studies, we included climatic parameters related to temperature and water availability, since these can be directly related to the physiological limitations of the study group. In addition, however, we also included relative elevation as a factor. This allows a direct comparison Abbreviations: absElev, absolute elevation; relElev, relative elevation; Bio5, max temperature of the warmest month; maxPET, maximum potential evapotranspiration; Bio14, precipitation of the driest month.
TA B L E 2 Summary of models regressing overall relative species richness. of mountain transects of different extensions and is somewhat independent of climate. We used it as an indicator of spatial processes, such as source-sink effects (Kessler, Hofmann, et al., 2011), overlaps between lowland and highland floras (Grytnes, 2003a), and island effects on small mountain ranges (Kessler & Kluge, 2022).

F I G U R E 3
The most important result of our study is that relative elevation is shown as the most powerful single predictor of species richness.
Previous studies of global elevational richness patterns of other groups of organisms, found climatic factors or absolute elevation to be dominant, although this cannot be directly compared with our study, since most of them did not calculate relative elevation.
However, in a study of fern richness using a similar number of transects as in the present study, and indeed, some of the same transects, , also found a predominance of unimodal richness patterns, but in contrast to the present study, the fern transects showed a very confused pattern with low explanatory power when aligned by relative elevation. Importantly, on short transects (<2000 m), richness either increased linearly (tropics) or decreased linearly (temperate regions) with increasing elevation, as also found by Khine et al. (2019). Instead of relative elevation, the temperature was found as the best explanatory factor, and thus,  concluded that fern richness is limited at low elevations by high temperatures and high elevations by low temperatures and that richness peaks under cool (15-18°C) and moist conditions as found at mid-elevations (1800-2000 m) along extensive tropical gradients. Considering that liverworts and ferns share a preference for cool and moist habitat conditions (Mandl et al., 2010), we here expected a similar result. Yet, this was not the case. In particular, liverwort richness showed hump-shaped patterns even on some mountains with limited elevational extent.
In addition to relative elevation, we also found that temperature played an important role in the models combining several explanatory factors. Interestingly, the majority of the individual gradients, had a hump-shaped relationship to Bio5, although the temperature of the richness peaks differed strongly. Temperature is a well-known factor driving species richness patterns both because it influences ecosystem productivity and hence resource availability for organisms (Chollet et al., 2014) and because many organisms have physiological limitations directly linked to high or low temperatures (e.g., Arris & Eagleson, 1989;Breitbarth et al., 2007;Reitzel et al., 2013). Several previous studies have found temperature to be an important predictor of liverwort species richness Henriques et al., 2016;Marline et al., 2020). In our study, we found maximum temperature of the warmest month (Bio5) to be the most important temperaturerelated variable, outranking both mean annual temperature (Bio1) and minimum temperature of the coldest month (Bio6). This indicates that liverwort richness appears to be more strongly limited by high than low temperatures, possibly because liverworts are poikilohydric, so that they cannot control evaporative water loss and depend on constant external water supply to maintain photosynthesis (Song, Zhang, et al., 2015). Thus, although liverworts are desiccation tolerant to varying degrees (Franks & Bergstrom, 2000;León-Vargas et al., 2006), active metabolism for growth inevitably requires available liquid water for turgescence of the cells (Proctor et al., 2007). In consequence, the growth and vitality of liverworts are reduced when they experience drought events (Song et al., 2012). Because water loss is closely linked to high temperatures, it is thus not surprising to find that these are related to liverwort richness.
The dependency of liverwort diversity on water availability is also reflected by the global trend of increasing species richness with increasing precipitation of the driest month (Bio14). For maxPET, which reflects the evaporative potential, we found a hump-shaped richness pattern, which may be interpreted as richness being limited by low energy availability and physiological limitations under cold conditions, and by high desiccation under high-temperature conditions.
So how can we explain the global mid-elevation richness peak in liverworts, as relElev performs best in explaining the elevational liverwort species richness pattern but cannot be interpreted as a directly acting driver of species richness? The maximum temperature of the warmest month (Bio5) in fact correlates with elevation, but a strict climate-based explanation as found for ferns cannot account for this pattern, since this would predict linear patterns on gradients of limited extent and a stronger predictive power of climatic factors in the analyses combining all elevational gradients. One might consider that liverwort richness is driven by climatic factors not sampled by us, such as the elevation of cloud condensation belts, which often occur at mid-elevation on extensive tropical gradients (Bruijnzeel et al., 2011;Hemp, 2006). However, it is unlikely that this is the missing factor, since condensation belts would occur at the top of gradients to a limited extent, and because their occurrence is correlated to temperature. Indeed, although they also did not have cloud cover as a factor,  still managed to explain fern richness by climatic factors.
The most likely explanation left is that liverwort richness is indeed somehow defined by factors or processes related to the relative position along the elevation gradient. Hump-shaped patterns of species richness have been postulated to occur when gradients are defined by strict barriers, in our case sea level and mountain top on most individual gradients (Colwell & Hurtt, 1994). In such a situation, species will overlap in the middle of the gradient, creating a humpshaped diversity pattern. This pattern was for a while explained via a null model, the mid-domain effect, which assumes no biological processes beyond the random placement of species (Colwell & Hurtt, 1994). However, this model has fundamental conceptual issues and has fallen out of favor (Currie & Kerr, 2008;Hawkins et al., 2005;Kluge et al., 2006), although it may influence patterns in subtle ways (Colwell et al., 2016).
Alternatively, there are a number of biologically realistic processes leading to hump-shaped elevational patterns. Largely, these are based on an accumulation of species at mid-elevations that arrive from low and high elevations (Grytnes, 2003a). This overlap zone can either be caused by non-self-sustaining sink populations (Kessler, Hofmann, et al., 2011) or by a more permanent overlap of different floras (Grytnes, 2003a). Since the ecological conditions of mountains change within short geographical distances, liverwort community composition changes considerably with elevation (Frahm & Gradstein, 1991;Gradstein et al., 2001). The mid-elevational richness peak may thus result from an overlap of species ranges from lower elevational zones with the ranges of species from high elevations, as previously found by Wolf (1993)  Finally, the models separating epiphytic and non-epiphytic liverwort species richness contained much less data (163 records from 16 and 90 records from nine locations, respectively) than the data set for the overall species richness (where published data often did allow the separation of epiphytic and non-epiphytic records), which probably hampered the detection of a clear global trend.
Furthermore, particularly non-epiphytic species richness of liverworts has been shown to be strongly associated with local-scale variables such as relative humidity, vegetation structure, and exposure (Mandl et al., 2009;Maul et al., 2020), which may obscure global trends especially with limited data. We thus do not further explore the results of these analyses.

| CON CLUS IONS
In this study, we for the first time compared worldwide data of elevational richness patterns of liverworts and found that a hump-shaped pattern with the highest richness at the middle of the elevational extent of each gradient is the predominant pattern. Compared to similar studies in other plant groups, this pattern appears to be more pronounced in liverworts than in ferns (Hernández-Rojas et al., 2021;Khine et al., 2019), woody angiosperms (Qian, 2018;Yue & Li, 2021), seed plants (Baniya et al., 2012;Zhang et al., 2016), herbs and trees (Vetaas et al., 2019), or vascular plants in general (Grytnes, 2003b). Whether this pattern truly arises from a great floristic turnover between high-altitude and lowland, remains to be confirmed by studies specifically assessing community structure. There is also a marked influence of climate-related parameters on liverwort richness, particularly via high temperatures, suggesting that liverwort richness in mountains is likely to be affected by temperature shifts associated with global warming.

ACK N OWLED G M ENTS
We thank the national authorities in Papua New Guinea for granting the permits to carry out fieldwork and collect voucher specimens. We

FU N D I N G I N FO R M ATI O N
Open Access funding enabled and organized by Projekt DEAL.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data set of this study (i.e. species numbers, elevational distribution, and (predicted) climate data as well as the original and determined GPS coordinates for each transect) is available from the DRYAD repository at https://doi.org/10.5061/dryad.4f4qr fjg6.